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Abstract: Cervus sichuanicus is a species of sika deer (Cervus nippon Group). To date, research has mainly focused on quantity 
surveying and behavior studies, with genetic information on this species currently deficient. To provide scientific evidence to assist in 
the protection of this species, we collected Sichuan sika deer fecal samples from the Sichuan Tiebu Nature Reserve (TNR) and 
extracted DNA from those samples. Microsatellite loci of bovine were used for PCR amplification. After GeneScan, the genotype 
data were used to analyze the genetic diversity and population structure of the Sichuan sika deer in TNR. Results showed that the 
average expected heterozygosity of the Sichuan sika deer population in TNR was 0.562, equivalent to the average expected 
heterozygosity of endangered animals, such as Procapra przewalskii. Furthermore, 8 of 9 microsatellite loci significantly deviated 
from the Hardy-Weinberg equilibrium and two groups existed within the Sichuan sika deer TNR population. This genetic structure 


may be caused by a group of Manchurian sika deer (Cervus hortulorum) released in TNR. 
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The Sichuan sika deer (Cervus sichuanicus) (Guo et 
al, 1978) is a species of the Cervus nippon group (Groves 
and Grubb, 2011). It is an endemic species in China and 
is listed as a class | endangered species, as well as an 
endangered species by the International Union for 
Conservation of Nature and Natural Resources (IUCN). 
Currently, the Sichuan sika deer is distributed in three 
unconnected areas located in northwest Sichuan Province 
and southwest Gansu Province. Compared to other spec- 
ies in China, Sichuan has the largest sika deer population 
in the wild, consisting of about 850 individuals (Guo, 
2000). 

In recent years, most research has focused on 
macroscopic aspects such as social behavior (Guo et al, 
1991), distribution and quantity surveying (Guo, 2000), 
activity rhythms (Liu et al, 2004), and food habits (Guo, 
2001). At the genetic level, Wu et al (2008) chose 16 
microsatellite loci to analyze the genetic diversity of Sika 
deer in China, including 24 samples of Sichuan sika deer. 
Lü et al (2006) and Liu et al (2003) studied the genetic 
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relationship of sika deer using mitochondrial DNA, 
which also had a small number (16) of Sichuan sika deer 
samples. These studies focused on the differences among 
species and genetic structure of sika deer in China. Due 
to the limited samples, however, the above data do not 
reflect genetic information in one species or population. 
Moreover, using traditional capture methods makes it 
difficult to obtain abundant samples to study genetic 
diversity and population structure within Sichuan sika 
deer populations. 

The Tiebu Nature Reserve (TNR) in the Zoige 
administrative district of northwest Sichuan has the 
largest sika deer population. According to monitoring 
data from TNR, there were about 1 000 Sichuan sika deer 
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in the reserve in 2012. TNR has been chosen as the main 
area to study Sichuan sika deer due to large population 
size and better protection situation (Li, 2010; Ning et al, 
2008; Qi et al, 2010). To date, however, no genetic 
information studies on TNR sika deer have been 
conducted and the genetic diversity and possible 
population substructures remain unknown. Previous 
studies have shown that the geographical distribution 
pattern of the population is related to geographical 
isolation (Boyce et al, 1999; Gavin et al, 1999; Scandura 
et al, 1998). However, it is unclear whether the rivers, 
roads, villages, and crop protection fences across the 
reserve affect gene flow or generate genetic structures 
within the 
hypothesized that geographic and artificial factors may 


sika deer population. Therefore, we 
affect the gene flow of sika deer, and thus established 
seven sampling areas in TNR to obtain sika deer fecal 
samples. Overall, we collected five hides, one tissue 
sample, and 149 fecal samples. After extracting DNA 
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from those samples, nine Bovinae microsatellite loci 
were amplified and the products were used for gene 
scanning. We then used the genotype data to analyze 
genetic diversity and population structure. The results of 
this study provide a scientific basis for the protection of 
the Sichuan sika deer. 


MATERIALS AND METHODS 


Study area and sampling collection 

Fecal samples were collected in seven areas in TNR 
(E102°46'-103°14’, N33°58'-34°16', about 206 922 km”), 
as shown in Figure 1. To avoid cross-contamination 
when sampling, we only collected one pellet from each 
fecal heap. We collected five hides, one tissue sample, 
and 149 fecal samples. The hide samples were collected 
from dead individuals in the wild. The only tissue sample 
came from an individual that had been killed by poachers. 
All samples were stored in 95% ethanol. 


Tiebu Nature Reserve 


Figure 1 Distribution of sampling areas in Tiebu Nature Reserve 


DNA extraction and amplification 

DNA was extracted using an E.Z.N.A HP Tissue 
DNA Maxi Kit and a Stool DNA kit (OMEGA) 
according to the manufacturer’s protocols. The former 
primers of nine Bovinae microsatellite loci (IDVGA29, 
BM4107, RT1, TGLA53, RM188, RT13, BM6506, 
BL42, CSSM019, ETH225, and CERVID14) were mo- 
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dified fluorescent groups. The PCR mixture contained 
1.5 uL genomic DNA (100 ng/pL), 1 uL of each pri- 
mer (10 pmol/uL), 0.5 uL BSA (20 mg/L), and 7.5 uL 
of Premix Ex Taq (Takara), with deionized water used to 
make the sample volume up to 15 uL. The PCR amp- 
lification was carried out using an initial denaturation at 
94 °C for 3 min, followed by 35 cycles at 94 °C for 30 s, 
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annealing for 30 s, primer extension at 72 °C for 45 s, 
and a final extension at 72 °C for 10 min. The 
products were preserved at 4 °C. One tissue sample 
was used as positive control and water as negative 


control to ensure that PCR was performed correctly 
and to avoid contamination. The PCR products were 
sent to Shanghai MAP Biotechnology Co., Ltd 
(Shanghai, China) for GeneScan. 


Table 1 Information on the nine microsatellite bovine loci used in this research 








f i A Annealing Sources of 
Locus Repeat motif Primer sequence Size range (bp) temperature (C) primers 

CCC ACA AGG TTA TCT ATC TCC AG 

IDVGA29 (AC)n 137-157 60 rN ea 
CCA AGA AGG TCC AAA GCA TCC AC (1996) 
AGC CCC TGC TAT TGT GTG AG f 

BM4107 (An (TC)n (TG)n 162-172 60 a 3 a 
ATA GGC TTT GCA TTG TTC AGG (1994) 
TGC CTT CTT TCA TCC AAC AA p l 

RT1 (GT)n 222-240 55 Wilson & White 
CAT CTT CCC ATC CTC TTT AC (1998) 
CAG CAG ACA GCT GCA AGA GTT AGC 

TGLAS3 (AC)n 176-216 50 Hoffmann et al 
CTT TCA GAA ATA GTT TGC ATT CAT GCA G (2004) 
GGG TTC ACA AAG AGC TGG AC 

RM188 (AC)n 141-151 50-52 paraa a 
GCA CTA TTG GGC TGG TGA TT (1994) 
GCC CAG TGT TAG GAA AGA AG : 

RTI (GT)n 293-324 60 ss 
CAT CCC AGA ACA GGA GTG AG (1994) 
GCA CGT GGT AAA GAG ATG GC 

BM6506 (AC)n 196-212 59-63 Moore etal 
AGC AAC TTG AGC ATG GCA C (1994) 
CAA GGT CAA GTC CAA ATG CC ; 

BL42 (AC)n 246-258 58-60 Bop al 
GCA TTT TTG TGT TAA TTT CAT GC (1994) 
TTG TCA GCA ACT TCT TGT ATC TTT ; : 

CSSMO019 (AC)n 138-156 52-58 Wilson & White 
TGT TTT AAG CCA CCC AAT TAT TTG (1998) 
GAT CAC CTT GCC ACT ATT TCCT 

ETH225 (AC)n 140-193 60 Steffen et al 
ACA TGA CAG CCA GCT GCT ACT (1993) 
TCT CTT GCG TCT CCT GCA TTG AC 

CERVID14 (AC)n 214-231 61-65 De Woody:etal 
AAT GGC ACC CAC TCC AGT ATT CTT C (1995) 


Species identification and population assessment 

We amplified all samples using two pairs of primers 
of the mitochondrial control region and compared the 
sequences in GenBank to ensure the samples in our 
research were from sika deer. We used GeneMapper v. 4 
(Applied Biosystems) to obtain the microsatellite 
fragment sizes and identified the heterozygotes and 
homozygotes. We did two repeats if the genotypes of the 
two repetitive PCRs were different; for example, if there 
was a heterozygote in one repeat but a homozygote in the 
other repeat. If there was a heterozygote in three of four 
repeats, we considered this site to be a heterozygous 
locus; the same for a homozygote locus. The probability 
of identity (PI) and the PI between siblings were 
calculated in Gimlet 1.3.3 (Valière, 2002; Liu & Yao, 
2013) using genotype data of the nine microsatellite loci. 
We used the following principle when using Mstools 
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(Park, 2001) to assess the population number: we 
allowed one genetic mismatch at one allele for one locus 
to be considered as two samples belonging to the same 
individual (Bellemain et al, 2005). To detect the validity 
of the above principle, we randomly selected 50 fecal 
samples (accounting for a third of the total number) and 
amplified the nine loci twice. At the same time, we 
amplified another two loci (ETH225 and CERVID14) 
of the selected 50 samples. We tested whether the 
influence of population number may be caused by the 
repeat amplification and increasing loci. 


Genetic diversity and population structure 

We computed the expected heterozygosity (Hg), 
observed heterozygosity (Ho), and polymorphism infor- 
mation content (PIC) using Mstools. Arlequin 3.1 (Exc- 
offier et al, 2005) and Fstat 2.9.3.2 (Goudet, 2001) were 
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used to test the Hardy-Weinberg equilibrium and Fis 
separately. The fixation index (Fsr) (Weir & Cockerham, 
1984) among sampling areas was tested in Arlequin 3.1, 
using 1000 random iterations to calculate confidence 
intervals (CJ). Software structure 2.3 (Pritchard et al, 
2000) was used to estimate the potential populations (K). 
The admixture model was used and run three times using 
10 replications of K ranging from 1 to 7. For those runs 
the following conditions were adopted: a burn-in period 
at 10 000 following 10 000 000 replicates of the MCMC. 
The AK calculated as per Evanno et al (2005) was used 
to estimate optimal K. The formulas were: 
AK=m|L(K+1)-2L(k)-L(K—1)|/s[L(K)]. We referred to 
the structure output InP(D) as L(K). m was the mean 
value and s was the standard deviation. Software 
CLUMPP Windows.1.1.2 (Jakobsson & Rosenberg, 2007) 
and destruct (Rosenberg, 2004) were used to visualize 
in the 


the individual coefficients of membership 


subpopulation calculated by structure 2.3. 


RESULTS 


Species identification and population assessment 

Comparing the obtained sequences of mitochondrial 
hypervariable region I from our research to the 
corresponding sequence in GenBank, all samples were 
obtained from sika deer. Thereinto, 143 of 155 samples 
belonged to the haplotypes of Sichuan sika deer and the 
remainder belonged to the haplotypes of Manchurian 
sika deer (Cervus hortulorum). 

Most of the Sichuan sika deer samples could be 
amplified at all the nine microsatellite loci, few samples 
(11) had a loss at one locus. Analyzing the united PI 
values for nine loci, we found that the discrimination 
power of nine loci was close to 100% (Figure 2). 
According to the rule of Bellemain, we assessed 76 sika 


deer individuals from 143 face samples. After the second 
amplification of the randomly selected samples, the 
population number did not change. The additional two 
loci (ETH225 and CERVID14) increased the population 
by one individual, accounting for 3% of the total number. 
The distribution of the 76 individuals in the samplings 
areas is shown in Table 3. 
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Figure 2 Multi-loci PI (biased and unbiased) and PZ (sibs) for 

sika deer in increasing order of single-locus values 
Genetic diversity and population structure 

In the Sichuan sika deer population in TNR, the 
mean number of alleles / locus was 6.56 and the Hg was 
0.562. The PIC values were higher than or close to 0.5, 
except for RT13 (0.129). All detected loci deviated from 
the Hardy-Weinberg equilibrium significantly (Figure 2) 
and the P value was 0. In addition to the locus RT13, the 
Fis values of the remaining eight loci were negative. The 
values of Hz ranged from 0.4292 to 0.6865 in the seven 
sampling areas and all Ho values were higher than those of 
Hy (Table 3). The Fsr among the seven sampling areas 
ranged from -—0.04781 to 0.05243. The highest Fr 
(0.05243) appeared between Area IV and VI but the lowest 
(-0.04781) appeared between Area I and III (Table 4). 


Table 2 Microsatellite diversity indices of Sichuan sika deer in Tiebu Nature Reserve 





Locus No. of alleles Hg Ho PIC Fis 
IDVGA29 5 0.587 0.722 0.516 —0.230 
BM4107 7 0.560 0.654 0.495 —0.169 
RTI 6 0.633 0.951 0.558 0.506 
CSSM019 6 0.635 0.926 0.560 0.462 
TGLAS3 8 0.643 0.899 0.571 0.401 
BL42 9 0.607 0.913 0.523 0.509 
RM188 6 0.656 0.738 0.607 0.125 
BM6506 8 0.604 0.974 0.519 —0.620 
RT13 4 0.134 0.026 0.129 0.810 
Mean 6.56 0.562 0.756 0.496 0.347 
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Table 3 Indices of genetic diversity of Sichuan sika deer in different sampling areas 





Sample area Sample size Loci typed Hg Ho 

I 10 9 0.4292 0.7160 
II 13 9 0.5499 0.8034 
Il 7 9 0.4831 0.7778 
IV 8 9 0.6804 0.7284 
V 20 9 0.6865 0.7444 
VI 11 9 0.4350 0.7063 
VI 1 9 0.5368 0.8571 


Table 4 F-statistics (Fsr, below diagonal) and significant differences among different areas 








Sample area I I M IV V VI VI 
I + = + + z + 
I —0.00493 - - — + - 
MI —0.04781 —0.03753 — — — — 
IV 0.03553 —0.0052 —0.00225 z + — 
V 0.05177 0.00789 0.01646 —0.01778 + 2 
VI —0.02823 —0.00142 —0.0334 0.03287 0.05243 + 
Vi 0.01389 —0.02825 —0.03837 —0.01869 0.01487 —0.00278 

+ denotes significant difference, — in contrast (95% confidence interval). 
Population structure 
Structure analysis showed that when K=2, the 4K DISCUSSION 


exhibited an obvious apex (Figure 3). The Fsr between 
subpopulations was 0.535. Individual coefficients of 
membership in subpopulations (Figure 4) showed that 
when K=2, the genetic amount of cluster 1 (Yellow) was 
larger and distributed in all sampling areas. In contrast, the 
genetic amount of cluster 2 (Blue) was less than cluster 1, 
and had a relatively smaller distribution area. Most of the 
blue genotype was distributed in Areas IV and V. 
80 


Figure 3 Value of AK as a function of K based on 10 runs 





I II I IV Vv VI VI 


Figure 4 Individual coefficients of membership in each 
sampling area 
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Population assessment 

In recent years, non-invasive sampling methods 
have been widely used in genetic studies of endangered 
animals (He et al, 2010; Wang et al, 2012). Because of 
the endangered status and fugacious feature of sika deer, 
non-invasive sampling methods (obtaining fecal samples) 
have become an important means of population 
assessment and individual identification. In this study, 
the stability of repeat amplification of nine microsatellite 
loci had no effect on assessing the population number, 
that is, our genotype data were reliable. After adding two 
loci, there was a 3% error between previous and post 
assessment results, which is much higher than 2% 
(Bellemain et al, 2005). Bellemain et al (2005) attributed 
such errors in population assessment to pseudogenes, 


pollution and missing data. 


Genetic diversity and differentiation 
Genetic diversity helps populations adapt to 
changing environments. With more variation, it is more 
likely that some individuals in a population will possess 
variations of alleles that are suited for the environment. 
The population will continue for more generations 
because of the success of these individuals (Khater et al, 


2011). In our research, the Hg of Sichuan sika deer 
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population was 0.562, close to the mean Hg (0.559) of 
sika deer in China, but much higher than the mean Hg for 
Sichuan sika deer (0.477) and slightly lower than the 
mean Hg for Manchurian sika deer (0.584) and south 
China sika deer (Jiangxi, 0.585; Zhejiang, 0.589) (Wu et 
al, 2008). Compared to other Ungulata, such as red deer 
(Cervus elaphus) (Hz=0.78, Hajji et al, 2007; Hp=0.804, 
Pérez-Espona et al, 2008; Hy=0.62—0.85, Zachos et al, 
2007), forest musk deer (Moschus _ berezovskii) 
(H,=0.8-0.9, Zhou et al, 2005), and Mongolian wild ass 
(Equus hemionus) (Hr=0.83, Kaczensky et al, 2011), 
Sichuan sika deer had a relatively lower Hg and was 
equivalent to Przewalski's Gazelle (Procapra przewalskii) 
(H:=0.552, Yang & Jiang, 2011) and milu (Elphurus 
davidianus) (Hz=0.46—0.54, Zeng et al, 2007). The Hgs 
values were different in each sampling area. Area I had 
the lowest Hg (0.4292), while Area V had the highest 
(0.6865). Polymorphism information content (PIC) is a 
measure of variation between microsatellite loci. It is 
generally considered that a PIC value greater than 0.5 is 
for highly polymorphic loci, a PIC value of less than 0.5 
and greater than 0.25 is for moderately polymorphic loci, 
and a PIC of less than 0.25 is for low polymorphic loci 
(Niu et al, 2001). In this study, the PIC values of six of 
the nine microsatellite loci were higher than 0.5, 
indicating highly polymorphic loci, while two were 
moderately polymorphic loci and one was a low 
polymorphic locus. Heterozygosity and PIC information 
indicated Sichuan sika deer in TNR had abundant genetic 
diversity. 

Hartl & Clark (1997) considered that population 
deviation from the Hardy Weinberg equilibrium was 
mainly due to small populations, not random mating, 
gene mutation, migration or other factors. The present 
study showed that nine microsatellite loci of the Sichuan 
sika deer population in TNR deviated significantly from 
the Hardy Weinberg equilibrium. The Sichuan sika deer 
population numbers in TNR grew steadily from 650 in 
2000 (Guo, 2000) to 1000 in 2012 (monitoring data from 
TNR). As a result, the deviation from the Hardy 
Weinberg equilibrium is unlikely attributed to a 
reduction in the population. The coefficient of genetic 
differentiation (Fsr), on the other hand, is an important 
indicator of genetic differentiation among subpopul- 
ations. According to Wright (1978), Fsrs values bet- 
ween 0-0.5 suggest no differentiation between sub- 
groups. The Fsrs values in our seven sampling areas were 
all less than 0.5, ranging from 0.04781 to 0.05243, and 
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showing negative value in some regions. This indicated 
that the genetic variation among the seven sampling 
areas mainly came from inside the sampling area, and the 
geographical barriers and human factors do not affect 
gene flow in Sichuan sika deer in TNR. In addition, due 
to the relatively closed geographical position of TNR, 
very few little sika deer were distributed in the periphery 
of the protected areas (Guo, 2000). Gene flow through 
immigration and emigration could not significantly affect 
the Sichuan sika deer populations to deviate from the 
Hardy Weinberg equilibrium. 

Microsatellite polymorphism studies of sika deer in 
China found that 16 microsatellite loci significantly 
deviated from the Hardy Weinberg equilibrium, but Hgs 
were higher than Hos (Wu et al, 2008). In the present 
study, Hos were significantly higher than the Hes. The 
Fis values of eight microsatellite loci were negative, 
showing heterozygosis excess. When a population 
experiences a reduction in its effective size, it generally 
develops heterozygosity excess at selectively neutral loci 
(Cormuet & Luikart, 1996). Xu et al (2009) used 
microsatellites to detect the population structure of 
Himalayan marmots and found heterozygosity excess 
due to historic population decline. Wu et al (2008) tested 
bottlenecks in four Chinese sika deer species and found 
that Sichuan sika deer had not experienced a population 
Therefore, the 
heterozygosity excess of Sichuan sika deer in TNR could 


bottleneck in recent history. 


not be caused by population reduction. 


Population structure 

A population composed of individuals from two 
different source populations will tend to have an excess 
of heterozygotes (Milkman, 1975). In testing population 
structure in the present study, we found that when K=2, 
AK exhibited a significant apex, indicating that two 
subpopulations existed within the Sichuan sika deer 
population. Considering that the Sichuan sika deer 
population had heterozygote excess and two subpo- 
pulations, we speculated that the sika deer in TNR might 
have two different sources populations. According to Wu 
(2004) and TNR staff, a captive-bred Manchurian sika 
deer group was bred and released in TNR in the late 
1960s. The mitochondrial data showed that haplotypes of 
Manchurian sika deer (unpublished data) do exist. 
Previous research has shown that introduced species 
often hybridize with local sibling species (Avise et al, 
1974; Gutiérrez, 1993; Perry & Goudet, 2001) and 
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panmictic admixtures can develop rapidly over a short 
period of time (Echelle & Connor, 1989; Echelle, 1987; 
1997). 
introduced sika deer have been shown to hybridize with 
native red deer (Senn & Pemberton, 2009). Therefore, 
there is a strong possibility that hybridization occurred 


Hybridization abounds in Cervidae, and 


between the two sika deer species in our research. The 
released Manchurian sika deer explained the existence of 
subpopulations in the Sichuan sika deer population in 
TNR. We speculate that the Manchurian sika deer 
genotype is the minority — the blue parts in the individual 
coefficients of membership, which were mainly 
distributed in Areas [V and V- because the number of 
introduced Manchurian sika deer is not clear and of the 
155 sika deer fecal samples examined, only 12 belonged 
to the haplotypes of Manchurian sika deer. The Sichuan 
sika deer population were in larger quantities than that of 


Manchurian sika deer. So, the yellow parts in Figure 4 
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